Conventional Molecular Dynamics Simulations Research Articles

Effect of quercetin on the protein-substrate interactions in SIRT6: Insight from MD simulations

SIRT6 is of interest for its promising effect in the treatment of aging-related diseases. Studies have shown quercetin (QUE) and its derivatives have varying degrees of effect on the catalytic effect of SIRT6. In the research, the effect of QUE on the protein-substrate interaction in the SIRT6-mediated mono-ADP ribosylation system was investigated by conventional molecular dynamics (MD) simulations combined with MM/PBSA binding free energy calculations. The results show that QUE can bind stably to SIRT6 with the binding energy of −22.8 kcal/mol and further affect the atomic interaction between SIRT6 and NAD+ (or H3K9), resulting in an increased affinity between SIRT6-NAD+ and decreased SIRT6-H3K9 binding capacity. At the same time, the binding of QUE can also alter some structural characteristics of the protein, with large shifts occurring in the residue regions involving the N-terminal (residues 1–27), Rossmann fold regions (residues 55–92), and ZBD (residues 164–179). Thus, QUE shows great potential as a scaffold for the design of novel potent SIRT6 modulators.

Understanding the Modulatory Role of E2012 on the γ-Secretase-Substrate Interaction.

Allosteric modulation plays a critical role in enzyme functionality and requires a deep understanding of the interactions between the active and allosteric sites. γ-Secretase (GS) is a key therapeutic target in the treatment of Alzheimer's disease (AD), through its role in the synthesis of amyloid β peptides that accumulate in AD patients. This study explores the structure and dynamic effects of GS modulation by E2012 binding, employing well-tempered metadynamics and conventional molecular dynamics simulations across three binding scenarios: (1) GS enzyme with and without L458 inhibitor, (2) the GS-substrate complex together with the modulator E2012 in two different binding modes, and (3) E2012 interacting with a C99 substrate fragment. Our findings reveal that the presence of L458 induces conformational changes that contribute to stabilization of the GS enzyme dynamics, previously reported as a key factor that allowed the resolution of the cryo-EM structure and the enhanced binding of E2012. Furthermore, we identified the most favorable binding site for E2012 within the GS-substrate complex, uncovering significant modulatory effects and a complex network of interactions that influence the position of the substrate for catalysis. In addition, we explore a potential substrate-modulator binding before the formation of the enzyme-substrate complex. The insights gained from our study emphasize the importance of these interactions in the development of potential therapeutic interventions that target the functionality of the GS enzyme in AD.

Generation of conformational ensembles of small molecules via surrogate model-assisted molecular dynamics

The accurate prediction of thermodynamic properties is crucial in various fields such as drug discovery and materials design. This task relies on sampling from the underlying Boltzmann distribution, which is challenging using conventional approaches such as simulations. In this work, we introduce surrogate model-assisted molecular dynamics (SMA-MD), a new procedure to sample the equilibrium ensemble of molecules. First, SMA-MD leverages deep generative models to enhance the sampling of slow degrees of freedom. Subsequently, the generated ensemble undergoes statistical reweighting, followed by short simulations. Our empirical results show that SMA-MD generates more diverse and lower energy ensembles than conventional MD simulations. Furthermore, we showcase the application of SMA-MD for the computation of thermodynamical properties by estimating implicit solvation free energies.

Unraveling the Ligand-Binding Sites of CYP3A4 by Molecular Dynamics Simulations with Solvent Probes.

Cytochrome P450 3A4 (CYP3A4) is one of the most important drug-metabolizing enzymes in the human body and is well known for its complicated, atypical kinetic characteristics. The existence of multiple ligand-binding sites in CYP3A4 has been widely recognized as being capable of interfering with the active pocket through allosteric effects. The identification of ligand-binding sites other than the canonical active site above the heme is especially important for understanding the atypical kinetic characteristics of CYP3A4 and the intriguing association between the ligand and the receptor. In this study, we first employed mixed-solvent molecular dynamics (MixMD) simulations coupled with the online computational predictive tools to explore potential ligand-binding sites in CYP3A4. The MixMD approach demonstrates better performance in dealing with the receptor flexibility compared with other computational tools. From the sites identified by MixMD, we then picked out multiple sites for further exploration using ensemble docking and conventional molecular dynamics (cMD) simulations. Our results indicate that three extra sites are suitable for ligand binding in CYP3A4, including one experimentally confirmed site and two novel sites.

Base pair dynamics, electrostatics, and thermodynamics at the LTR-III quadruplex:duplex junction

G-quadruplexes (GQs) play key regulatory roles within the human genome and have also been identified to play similar roles in other eukaryotes, bacteria, archaea, and viruses. Human immunodeficiency virus 1, the etiological agent of acquired immunodeficiency syndrome, can form two GQs in its long terminal repeat (LTR) promoter region, each of which act to regulate viral gene expression in opposing manners. The major LTR GQ, called LTR-III, is a distinct hybrid GQ containing a 12-nucleotide duplex loop attached to the quadruplex motif. The resulting quadruplex:duplex junction (QDJ) has been hypothesized to serve as a selective drug targeting site. To better understand the dynamics of this QDJ, we performed conventional and enhanced-sampling molecular dynamics simulations using the Drude-2017 force field. We observed unbiased and reversible formation of additional base pairs in the QDJ, between Ade4:Thy14 and Gua3:Thy14. Both base pairs were electrostatically favored, but geometric constraints within the junction may drive the formation of, and preference for, the Ade4:Thy14 base pair. Finally, we demonstrated that the base pairs are separated only by small energy barriers that may enable transitions between both base-paired states. Together, these simulations provide new insights into the dynamics, electrostatics, and thermodynamics of the LTR-III QDJ.

MD investigation on the differences in the dynamic interactions between the specific ligand azamulin and two CYP3A isoforms, 3A4 and 3A5

The unmarked potential drug molecule azamulin has been found to be a specific inhibitor of CYP3A4 and CYP3A5 in recent years, but this molecule also shows different binding ability and affinity to the two CYP3A isoforms. In order to explore the microscopic mechanism, conventional molecular dynamics (MD) simulation methods were performed to study the dynamic interactions between two isoforms and azamulin. The simulation results show that the binding of the ligand leads to different structural properties of two CYP3A proteins. First of all, compared with apo-CYP3A4, the binding of the ligand azamulin can lead to changes in the structural flexibility of CYP3A4, i.e., holo-CYP3A4 is more flexible than apo-CYP3A4. The structural changes of CYP3A5 are just the opposite. The ligand binding increases the rigidity of CYP3A5. Furthermore, the representative structures of the production phase in the MD simulation were in details analyzed to obtain the microscopic interactions between the ligand azamulin and two CYP3A isoforms at the atomic level. It is speculated that the difference of composition and interaction of the active sites is the fundamental cause of the change of structural properties of the two proteins. Communicated by Ramaswamy H. Sarma

Role of Protonation States in the Stability of Molecular Dynamics Simulations of High-Resolution Membrane Protein Structures.

Classical molecular dynamics (MD) simulations provide unmatched spatial and time resolution of protein structure and function. However, the accuracy of MD simulations often depends on the quality of force field parameters and the time scale of sampling. Another limitation of conventional MD simulations is that the protonation states of titratable amino acid residues remain fixed during simulations, even though protonation state changes coupled to conformational dynamics are central to protein function. Due to the uncertainty in selecting protonation states, classical MD simulations are sometimes performed with all amino acids modeled in their standard charged states at pH 7. Here, we performed and analyzed classical MD simulations on high-resolution cryo-EM structures of two large membrane proteins that transfer protons by catalyzing protonation/deprotonation reactions. In simulations performed with titratable amino acids modeled in their standard protonation (charged) states, the structure diverges far from its starting conformation. In comparison, MD simulations performed with predetermined protonation states of amino acid residues reproduce the structural conformation, protein hydration, and protein-water and protein-protein interactions of the structure much better. The results support the notion that it is crucial to perform basic protonation state calculations, especially on structures where protonation changes play an important functional role, prior to the launch of any conventional MD simulations. Furthermore, the combined approach of fast protonation state prediction and MD simulations can provide valuable information about the charge states of amino acids in the cryo-EM sample. Even though accurate prediction of protonation states in proteinaceous environments currently remains a challenge, we introduce an approach of combining pKa prediction with cryo-EM density map analysis that helps in improving not only the protonation state predictions but also the atomic modeling of density data.

The Role of Force Fields and Water Models in Protein Folding and Unfolding Dynamics.

Protein folding is a fascinating, not fully understood phenomenon in biology. Molecular dynamics (MD) simulations are an invaluable tool to study conformational changes in atomistic detail, including folding and unfolding processes of proteins. However, the accuracy of the conformational ensembles derived from MD simulations inevitably relies on the quality of the underlying force field in combination with the respective water model. Here, we investigate protein folding, unfolding, and misfolding of fast-folding proteins by examining different force fields with their recommended water models, i.e., ff14SB with the TIP3P model and ff19SB with the OPC model. To this end, we generated long conventional MD simulations highlighting the perks and pitfalls of these setups. Using Markov state models, we defined kinetically independent conformational substates and emphasized their distinct characteristics, as well as their corresponding state probabilities. Surprisingly, we found substantial differences in thermodynamics and kinetics of protein folding, depending on the combination of the protein force field and water model, originating primarily from the different water models. These results emphasize the importance of carefully choosing the force field and the respective water model as they determine the accuracy of the observed dynamics of folding events. Thus, the findings support the hypothesis that the water model is at least equally important as the force field and hence needs to be considered in future studies investigating protein dynamics and folding in all areas of biophysics.

Implementation and Validation of an OpenMM Plugin for the Deep Potential Representation of Potential Energy.

Machine learning potentials, particularly the deep potential (DP) model, have revolutionized molecular dynamics (MD) simulations, striking a balance between accuracy and computational efficiency. To facilitate the DP model's integration with the popular MD engine OpenMM, we have developed a versatile OpenMM plugin. This plugin supports a range of applications, from conventional MD simulations to alchemical free energy calculations and hybrid DP/MM simulations. Our extensive validation tests encompassed energy conservation in microcanonical ensemble simulations, fidelity in canonical ensemble generation, and the evaluation of the structural, transport, and thermodynamic properties of bulk water. The introduction of this plugin is expected to significantly expand the application scope of DP models within the MD simulation community, representing a major advancement in the field.

Unraveling the allosteric mechanisms of prolyl endopeptidases for celiac disease therapy: Insights from molecular dynamics simulations

Prolyl endopeptidases (PEP) from Sphingomonas capsulata (sc) and Myxococcus xanthus (mx) selectively degrade gluten peptides in vitro, offering a potential therapeutic strategy for celiac disease. However, the mechanisms governing the interaction of these enzymes with their substrates remain unclear. In this study, conventional molecular dynamics simulations with a microsecond timescale and targeted molecular dynamics simulations were performed to investigate the native states of mxPEP and scPEP enzymes, as well as their allosteric binding with a representative substrate, namely, Z-Ala-Pro-p-nitroanilide (pNA). The simulations reveal that the native scPEP is in an open state, while the native mxPEP is in a closed state. When pNA approaches a closed mxPEP, it binds to an allosteric pocket located at the first and second β-sheet of the β-propeller domain, inducing the opening of this enzyme. Neither enzyme is active in the open or partly-open states. Enzymatic activity is enabled only when the catalytic pocket in the closed state fully accommodates the substrates. The internal capacity of the catalytic pocket of PEP in the closed state determines the maximum size of the gluten peptides that the enzymes can catalyze. The present work provides essential molecular dynamics information for the redesign or engineering of PEP enzymes.

Modeling and Simulation of the NMDA Receptor at Coarse-Grained and Atomistic Levels.

N-Methyl-D-aspartate (NMDA) receptors are glutamate-gated excitatory channels that play essential roles in brain functions. While high-resolution structures were solved for an allosterically inhibited form of functional NMDA receptor, other key functional states (particularly the active open-channel state) have not yet been resolved at atomic resolutions. To decrypt the molecular mechanism of the NMDA receptor activation, structural modeling and simulation are instrumental in providing detailed information about the dynamics and energetics of the receptor in various functional states. In this chapter, we describe coarse-grained modeling of the NMDA receptor using an elastic network model and related modeling/analysis tools (e.g., normal mode analysis, flexibility and hotspot analysis, cryo-EM flexible fitting, and transition pathway modeling) based on available structures. Additionally, we show how to build an atomistic model of the active-state receptor with targeted molecular dynamics (MD) simulation and explore its energetics and dynamics with conventional MD simulation. Taken together, these modeling and simulation can offer rich structural and dynamic information which will guide experimental studies of the activation of this key receptor.

Conformational Selectivity of ITK Inhibitors: Insights from Molecular Dynamics Simulations.

Interleukin-2-inducible T-cell kinase (ITK) regulates the response to T-cell receptor signaling and is a drug target for inflammatory and immunological diseases. Molecules that bind preferentially to the active form of ITK have low selectivity between kinases, whereas those that bind preferentially to the inactive form have high selectivity for ITK. Therefore, computational methods to predict the conformational selectivity of compounds are required to design highly selective ITK inhibitors. In this study, we performed absolute binding free-energy perturbation (ABFEP) simulations for 11 compounds on both active and inactive forms of ITK, and the calculated binding free energies were compared with experimental data. The conformational selectivity of 10 of the 11 compounds was correctly predicted using ABFEP. To investigate the mechanism underlying the stabilization of the active and inactive structures by the compounds, we performed extensive, conventional molecular dynamics simulations, which revealed that the compound-induced stabilization of the P-loop and linkage of conformational changes in L489, V419, F501, and M410 upon compound binding were critical factors. A guideline for designing inactive-form binders is proposed based on these key structural factors. The ABFEP and the created guidelines are expected to facilitate the discovery of highly selective ITK inhibitors.

Uncovering the Structural and Binding Insights of Dual Inhibitors Simultaneously Targeting Two Distinct Sites on EGFR Kinase.

Epidermal growth factor receptor (EGFR) is the first growth factor receptor identified in normal cells that is related to the receptor tyrosine kinase, which causes regular cell division. A point mutation in EGFR intracellular kinase domain forces the abnormal cell divisions throughout time, leading to non-small cell lung cancer (NSCLC) transformation. Thus, competitive inhibitors that bind to the ATP binding pocket have been developed as a targeted therapy for NSCLC. The third-generation kinase inhibitor Osimertinib is currently playing a very vital role in the treatment of NSCLC. However, it is not effective toward the C797S kinase domain mutation. For this reason, fourth-generation kinase noncompetitive inhibitors are introduced which work through binding to an allosteric pocket near the ATP binding region and act as a better binding agent for this mutated kinase domain. However, the problem is that these single fourth-generation kinase inhibitors may not be as effective as a single agent. The aim of this work was to apply combinations of these two inhibitors together in different binding regions of EGFR without overlapping the resistance mechanism to obtain the key direct and indirect interactions occurring between them. Moreover, the free energy of dissociation of an inhibitor from its binding sites in the presence of a second inhibitor immobilized in another binding site was also the focus of the study. To realize this aim, we performed conventional molecular dynamics simulations and principal component analysis and dynamic cross-correlation matrices along with umbrella sampling. Our results demonstrated that binding of dual inhibitors triggered conformational changes of the protein more toward the inactive state. Furthermore, allosteric inhibitors bound more strongly to protein kinase EGFR than the orthosteric inhibitors in the presence of dual inhibitors. Finally, the binding mechanism and important hydrogen-bonding residues during unbinding of the inhibitors were fully elucidated. This study provides insight into the binding of the receptor-orthosteric inhibitor-allosteric inhibitor, which can be helpful for further design of novel inhibitors that have a better inhibitory action.

Molecular mechanism of calcitriol enhances membrane water permeability

Helicobacter pylori (H. pylori) exhibits a unique membrane lipid composition, including dimyristoyl phosphatidylethanolamine (DMPE) and cholesterol, unlike other Gram-negative bacteria. Calcitriol has antimicrobial activity against H. pylori, but cholesterol enhances antibiotics resistance in H. pylori. This study explored the changes in membrane structure and the molecular mechanisms of cholesterol/calcitriol translocation using well-tempered metadynamics (WT-MetaD) simulations and microsecond conventional molecular dynamics (CMD) simulations. Calcitriol facilitated water transport across the membrane, while cholesterol had the opposite effect. The differing effects might result from the tail 25-hydroxyl group and a wider range of orientations of calcitriol in the DMPE/dimyristoyl phosphatidylglycerol (DMPG) (3:1) membrane. Calcitriol moves across the bilayer center without changing its orientation along the membrane Z-axis, becomes parallel to the membrane surface at the membrane-water interface, and then rotates approximately 90° in this interface. The translocation mechanism of calcitriol is quite different from the flip-flop of cholesterol. Moreover, calcitriol crossed from one layer to another more easily than cholesterol, causing successive perturbations to the hydrophobic core and increasing water permeation. These results improve our understanding of the relationship between cholesterol/calcitriol concentrations and the lipid bilayer structure and the role of lipid composition in water permeation.

A novel constant pressure control method for nucleation boiling of CO2/lubricant in molecular dynamics simulation

Conventional molecular dynamics simulations of nucleation boiling typically employ the fixed volume method, resulting in a non-constant pressure during the boiling process. In this paper, we propose an additional pressure control method, placing a cover plate at the top of the liquid phase to achieve liquid phase pressure control during boiling. The novel method yields a pressure control error within ±5 %. CO2 and CO2/PEC4 lubricating oil mixture are used to investigate the differences in heat transfer characteristics between non-constant and constant pressure nucleation boiling. Furthermore, the effect of PEC4 lubricating oil on CO2 nucleation boiling under a fixed pressure condition is explored. The results show that both the bubble nucleation and growth time of CO2 and the CO2/PEC4 mixture are prolonged under constant pressure. This causes a compression in CO2 molecular spacing, thereby reducing CO2 diffusion and mass transfer capabilities. Heat transfer from the solid to the liquid, and within the liquid itself, is diminished, while the distribution of PEC4 molecules within the liquid phase becomes more even. The migration rate of PEC4 molecules in the liquid phase slows down. Under the constant pressure model, PEC4 increases the mass transfer resistance of CO2. As the lubricant mass fraction rises, The nucleation and growth rate of CO2 bubbles decline. In comparison to CO2, PEC4 molecules exhibit a stronger interaction energy with the Cu wall surface. During boiling, some PEC4 molecules accumulate on the Cu surface, forming an oil film that impedes the direct contact between CO2 molecules and the Cu wall.

GaMD simulations as an alternative in the TFE-water mixture description

Context: 2,2,2-Trifluoroethanol has been widely used to study the structure and dynamic properties of intrinsically disordered proteins. Experimentally, it is known that TFE-water mixtures stabilize secondary structures of IDPs, and therefore, it allows the studying of conformational ensembles of these proteins. In the last decades, molecular dynamic simulations have helped study the IDPs’ conformational ensemble. Unfortunately, conventional MD requires very long simulation times to describe the properties of IDPs. Therefore, a variety of accelerated sampling techniques have been developed and employed. The TFE-water mixture arrangement description through MD has faced substantial difficulties since emulating the TFE nanocrowding at certain TFE:H_2O ratios (around 15–40% of TFE). In this work, we determine the most suitable conditions that reproduce experimentally reported properties of TFE-water mixtures. We compared the employment of conventional MD and GaMD simulations and various water parameters. Our results show that the combination of parameters that better reproduce the experimental information is the combination of the TIP4PD water model and GaMD simulations. Therefore, these conditions help accurately describe the structural ensemble of IDPs in TFE-water mixtures. Methods: Conventional MD and GaMD simulations were performed under AMBER 18 software. The TFE and water molecules were described using GAFF2 and a variety of water models, such as TIP3P, TIP4P2005, TIP4PD, and TIP5P, respectively. The systems were simulated a 100 ns at 298 K.Graphical abstract

Activation Loop Plasticity and Active Site Coupling in the MAP Kinase, ERK2

Previous studies of the protein kinase, ERK2, using NMR and hydrogen-exchange measurements have shown changes in dynamics accompanying its activation by phosphorylation. However, knowledge about the conformational motions involved is incomplete. Here, we examined ERK2 using long conventional molecular dynamics (MD) simulations starting from crystal structures of phosphorylated (2P) and unphosphorylated (0P) forms. Individual trajectories were run for (5 to 25) μs, totaling 727 μs. The results show unexpected flexibility of the A-loop, with multiple long-lived (>5 μs) conformational states in both 2P- and 0P-ERK2. Differential contact network and principal component analyses reveal coupling between the A-loop fold and active site dynamics, with evidence for conformational selection in the kinase core of 2P-ERK2 but not 0P-ERK2. Simulations of 2P-ERK2 show A-loop states corresponding to restrained dynamics within the N-lobe, including regions around catalytic residues. One A-loop conformer forms lasting interactions with the L16 segment, leading to reduced RMSF and greater compaction in the active site. By contrast, simulations of 0P-ERK2 reveal excursions of A-loop residues away from the C-lobe, leading to greater active site mobility. Thus, the A-loop in ERK2 switches between distinct conformations that reflect coupling with the active site, possibly via the L16 segment. Crystal packing interactions suggest that lattice contacts with the A-loop may restrain its structural variation in X-ray structures of ERK2. The novel conformational states identified by MD expand our understanding of ERK2 regulation, by linking the activated state of the kinase to reduced dynamics and greater compaction surrounding the catalytic site.

Unraveling pyrolysis mechanisms of lignin dimer model compounds: Neural network-based molecular dynamics simulation investigations

Unraveling the pyrolysis mechanisms of lignin is of great importance for effective lignin utilization. However, both conventional quantum mechanical method and molecular dynamics (MD) simulations based on empirical potentials, e.g., the ReaxFF, cannot guarantee simultaneous accuracy and efficiency in exploring lignin pyrolysis mechanisms. Here, we developed a neural network (NN) force field for four types of representative lignin dimers based on 56,164 structures at the level of density functional theory (DFT). Such NN model was demonstrated to possess more accurate predictive power in both energetic and atomistic force properties than ReaxFF and was far more efficient than QM method. The MD simulations with the NN model showed that entire pyrolysis process can be divided into three stages dependent on the reaction temperature: homolytic cleavage of C–C bond and C-O bond, the changes from methoxy group to hydroxy, methyl, and formyl group, ring-opening reactions and the ring structure evolutions, and some novel reactions. The NN-based MD has provided us an atomic-level understanding of lignin pyrolysis and enabled us to build a comprehensive lignin pyrolysis network by means of the reactions disclosed based on the NN-based MD, which provided a useful guide for taking advantage of lignin more effectively.

Unraveling Tribochemistry and Self-Lubrication Mechanism of Polytetrafluoroethylene by Reactive Coarse-Grained Molecular Dynamics Simulations.

Lubrication of polymeric materials generally involves processes of atomic-scale chemical bond forming/breaking at the interface and mesoscale chain reorientation, disentanglement, and so forth. However, it is difficult to describe the important aspects of tribochemical reactions by conventional coarse-grained molecular dynamics (CGMD) simulations. Here, reactive CGMD simulations were conducted based on the ReaxFF force field to study the tribochemical interactions between polytetrafluoroethylene (PTFE) and iron. The chemical bond forming/breaking between the molecular chain and countersurface was fitted through the bond dissociation energies of specific reaction sites from all-atom ReaxFF-MD simulations. This enabled a quantitative description of tribochemical reactions in a macromolecule system. First, the number of anchoring bonds between PTFE molecules and the countersurface showed a strong correlation with the friction coefficient. The shearing process induced breaking of the interfacial anchoring bonds as well as chain disentanglement in the matrix, which consequently led to ordering reorientation of molecular chains toward sliding direction and hence decrease of friction. Second, two competitive factors were clarified to affect polymer friction with varying temperatures. The decrease of interfacial anchoring reactivity and molecular chain mobility at low temperature prohibited reorientation of molecular chains and increased the friction coefficient. On the other hand, the hardening of PTFE and the reduction in effective contact area at low temperatures decreased the friction coefficient. This led to a turning point with a maximum friction coefficient around 100 K. These results shed light on the essential role of tribochemical reactions on polymer lubrication, especially under low temperatures, which provides design guidance of polymeric lubrication systems for engineering applications.

Accelerating Molecular Dynamics Simulations for Drug Discovery.

Accurate prediction of ligand binding thermodynamics and kinetics is crucial in drug design. However, it remains challenging for conventional molecular dynamics (MD) simulations due to sampling issues. Gaussian accelerated MD (GaMD) is an enhanced sampling method that adds a harmonic boost to overcome energy barriers, which has demonstrated significant benefits in exploring protein-ligand interactions. Especially, the ligand GaMD (LiGaMD) applies a selective boost potential to the ligand nonbonded potential energy, significantly improving sampling for ligand binding and dissociation. Furthermore, a selective boost potential is applied to the potential of both ligand and protein residues around binding pocket in LiGaMD2 to further increase the sampling of protein-ligand interaction. LiGaMD and LiGaMD2 simulations could capture repetitive ligand binding and unbinding events within microsecond simulations, allowing to simultaneously characterize ligand binding thermodynamics and kinetics, which is expected to greatly facilitate drug design. In this chapter, we provide a brief review of the status of LiGaMD in drug discovery and outline its usage.